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ABSTRACT 


Short  period  P  waves  of  four  presumed  Soviet  explosions  from  Eastern  Kazakh 
are  examined  at  five  teleseismic  arrays:  LASA,  YKA,  OONW,  WRA  and  GBA.  Trans¬ 
fer  functions  to  shape  the  lower  magnitude  to  the  highest  magnitude  event  were  com¬ 
puted  at  each  array  to  eliminate  transmission  path  effects  from  source  to  receiver. 

The  transfer  functions  are  locally  consistent  from  sensor  to  sensor  at  each  array,  but 
show  considerable  global  variations  from  array  to  array.  This  suggests  that  there  are 
strong  azimuthal  variations  in  the  source  radiation  of  the  events,  due  to  complex 
scattering  of  the  signals  at  the  test  site.  At  LASA,  the  observed  transfer  functions 
can  be  explained  by  using  explosion  source  functions  of  Blake  and  Haskell.  The  assumed 
source  parameters  are  scaled  functions  of  the  yield  of  each  event,  which  is  estimated 
from  an  empirical  amplitude  yield  relationship. 


Accepted  for  the  Air  Force 

Joseph  R,  Waterman,  Lt.  Col.,  USAF 

Chief,  Lincoln  Laboratory  Project  Office 
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I.  INTRODUCTION 


In  a  previous  Technical  Note,  Filson^  discussed  the  short  period  spectra 
from  five  presumed  nuclear  explosions  from  Eastern  Kazakh  recorded  at  five  world 
wide  arrays.  His  analysis  consisted  in  fitting  the  corrected  displacement  spectrum  of 
each  event  at  each  array  to  a  scaled  explosion  source  model  given  by  Haskell.  In 
order  to  use  the  observed  data,  corrections  had  to  be  made  for  instrument  response, 
crustal  and  mantle  layering,  and  average  attenuation  effects  in  the  Earth.  In  spite  of 
the  uncertainties  in  estimating  these  quantities,  Filson  was  able  to  demonstrate  a 
systematic  shift  of  the  displacement  frequency  content  towards  lower  frequencies  with 
increasing  magnitude  as  predicted  by  Haskell's  explosive  source  model. 

In  this  report  a  different  technique  is  used  to  verify  the  apparent  seismic 
scaling  effect  with  magnitude  found  by  Filson.  Using  four  of  the  events  analyzed  by 
Filson,  transfer  functions  are  computed  in  time  which  shape  each  of  the  lower  magni¬ 
tude  events  to  the  largest  event  at  each  array.  For  each  pair  of  events  this  is  equiva¬ 
lent  to  taking  the  spectral  ratio  of  the  largest  P  wave  over  a  smaller  P  wave,  both 
recorded  at  the  same  site.  Such  transfer  functions  cancel  out  all  the  unknown  trans¬ 
mission  path  effects  and  the  instrument  response  at  each  site,  and  therefore  can  be 
interpreted  as  spectral  ratios  of  the  source  radiations  of  the  events. 

In  computing  transfer  functions  at  a  given  array,  one  can  use  either  the  data 
from  individual  sensors  or  steered  beams  for  the  entire  array.  Whereas  steered  beams 
generally  have  better  signal  to  noise  ratios  than  signal  sensors,  the  self-consistency  of 
the  data  can  be  tested  by  computing  transfer  functions  at  each  sensor.  If  such  functions 
are  coherent  across  the  entire  array,  then  they  can  be  reliably  interpreted  in  terms  of 
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source  functions.  The  first  part  of  the  data  analysis  shows  that  quite  consistent  trans¬ 
fer  functions  are  obtained  in  the  individual  sites  of  each  array,  in  spite  of  large  fluc¬ 
tuations  of  signal  shape  and  amplitude.  Thus  no  apparent  differences  in  radiation 
patterns  can  be  detected  within  each  array. 

Array  beams  of  each  event  are  shown  in  the  second  part  of  the  data  analysis. 
Transfer  functions  are  computed  from  this  data  at  each  array.  For  reasons  described 
in  that  section  the  LASA  transfer  functions  are  considered  most  reliable  for  interpre¬ 
tation  in  terms  of  source  functions. 

Two  sets  of  source  functions  are  calculated  for  the  events,  using  models 
given  by  Haskell  and  Blake.  Yields  for  the  four  events  are  estimated  from  an  empirical 
magnitude -yield  formula  for  explosions  in  hard  rock.  From  the  yields  the  Haskell  and 
Blake  solutions  are  obtained  by  scaling  the  yield  dependent  parameters. 

Finally,  theoretical  transfer  functions  are  computed  from  the  two  sets  of 
source  functions,  and  they  are  compared  to  the  LASA  transfer  functions. 
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II.  DESCRIPTION  OF  DATA  AND  ARRAYS 


The  short  period  data  used  in  this  analysis  was  recorded  at  five  arrays 
azimuthally  distributed  about  the  Soviet  test  site  in  Eastern  Kazakh.  These  arrays 
are  the  following:  (LASA),  the  Large  Aperture  Seismic  Array  in  Montana,  USA;  (GBA), 
near  Gauribidanur,  India;  (WRA),  near  Warramunga,  Australia;  (YKA),  near  Yellow¬ 
knife,  Canada;  and  (OONW)  near  Oslo,  Norway.  In  Figure  1  the  location  of  each  array 
is  shown  relative  to  the  Soviet  test  site  denoted  here  by  (STS).  The  epicentral  distance 
from  the  test  site  to  each  array  varies  from  ~36°  at  OONW  to  ~89°at  WRA. 

The  geometries  of  the  arrays  are  shown  in  Figures  2  and  3.  Clearly  LASA , 
with  an  aperture  of  200  km.  is  an  order  of  magnitude  larger  than  the  other  arrays. 

GBA,  YKA  and  WRA  are  arrays  operated  by  the  United  Kingdom,  and  OONW  was  a 
temporary  array  set  up  by  the  United  States,  but  now  not  operating.  These  smaller 
arrays  each  consist  of  two  orthogonal  arms  of  evenly  spaced  seismometers,  the  arms 
being  20  km.  long  or  less. 

The  short  period  displacement  response  of  the  seismometers  contained  in  the 
arrays  is  shown  in  Figure  4.  Converting  these  two  curves  to  velocity  response  one 
obtains  a  curve  for  LASA  and  UK  whi  ch  is  flat  from  about  1  to  4  Hz  and  an  OONW 
curve  which  is  flat  from  about  1  to  2  Hz  tapering  off  slowly  above  2  Hz. 

These  first  four  figures  are  taken  from  Filson's  report. 1 

The  data  consists  of  the  short  period  teleseismic  P  waves  from  four  presumed 
explosions  from  Eastern  Kazakh  recorded  at  the  arrays  described  above.  Table  1  is  a 
list  of  the  four  events  arranged  in  order  of  increasing  recorded  amplitude  at  most  of 
the  five  arrays.  The  event  numbers  refer  to  the  events  studied  by  Filson^.  The  data 
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Fig.  1.  World  map  showing  the  array  locations  relative  to  the  Soviet  test  site  (STS). 
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Fig.  2.  LAS  A  geometry.  Center  seismometer  is  at  46°  41 '19.0”  N,  106o13'20.0”  W. 


5 


62*35'N 

62*  30'  N 

1  1 

5km 

YELLOWKNIFE  ARRAY 

CANADA 

62*25'N 

1 

1 

1 14°  45 W 

114*  30W 

Fig.  3.  Geometry  of  the  United  Kingdom  arrays  and  OONW. 
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Fig.  4.  Shapes  of  the  response  curves  of  the  short  period 
vertical  instruments  of  the  arrays. 
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tabulated  without  asterisks  for  Events  1,  2  and  5  were  obtained  from  the  Earthquake 

Data  Report  published  by  the  USCGS.  Event  4,  however,  was  not  reported  by  the 

USCGS.  In  order  to  locate  this  event  relative  to  the  other  three,  first  motions  were 

picked  by  the  author  at  the  LASA  subarrays  for  all  four  events.  For  each  event,  a 

least  squares  plane  wave  was  calculated  for  the  subarray  arrival  times  using  the 
2 

Analysis  Console  at  Lincoln  Laboratory.  The  locations  obtained  are  listed  in  the 
last  column  of  Table  1. 


Event 

Date 

Origin  Time 

Location 

Location 

(LASA) 

h 

m 

s 

1 

5.3 

Sept  22,  *67 

5 

03 

58 

50.  ON 

77.  6E 

51.  3N 

79.  9  E* 

2 

5.3 

Sept  16,  '67 

4 

03 

58 

50.  ON 

77.  8  E 

50.  9N 

78.  4E* 

4 

— 

Aug  4,  '67* 

6 

57 

* 

00 

m 

50.  8N 

78.  IE* 

5 

5.7 

Oct  17,  *67 

5 

03 

58 

49.  8N 

78.  IE 

50.  3N 

78.  4E* 

Presumed  Explosions  from  E.  Kazakh.  Data  with  asterisks (*)  were  estimated  from 
LASA  data  by  the  author.  Data  without  asterisks  were  obtained  from  the  USCGS. 

Figure  5  shows  the  USCGS  and  LASA  locations  of  the  events  in  Eastern 
Kazakh.  As  expected  the  LASA  locations  are  not  as  tightly  grouped  as  the  USCGS  loca¬ 
tions.  However,  Event  4  appears  to  be  located  close  to  the  other  three  events.  For 
the  purposes  of  this  paper  all  events  are  assumed  to  have  nearly  identical  epicenters. 
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Fig.  5.  Map  of  the  Soviet  test  site  region  in  Eastern  Kazakh. 
Open  circles  are  USCGS  locations  and  solid  circles  are  first 
motion  LASA  locations  for  the  events  studied. 
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III.  CONSISTENCY  OF  SOURCE  RADIATION  TO  EACH  ARRAY 


A  crucial  question  is  whether  or  not  the  observed  seismic  radiation  is 
locally  consistent  over  the  aperture  of  each  array.  Figures  6  and  7  show  seismograms 
from  the  largest  and  smallest  arrays  of  this  study,  i.  e.  LASA  and  OONW.  At  LASA 
the  subarray  sum  traces  are  shown  and  at  OONW  the  individual  seismometer  traces 
are  displayed.  In  each  figure  one  can  see  considerable  variation  of  signal  shape  and 
duration.  These  variations  are  not  random  however,  but  are  quite  repeatable  from 
event  to  event  at  each  site.  This  suggests  that  the  teleseismic  radiation  arriving  at 
each  seismometer  is  severely  altered  by  scattering  in  the  crust  and  upper  mantle  under 
each  recording  site.  At  OONW  Events  1  and  2  have  identical  shapes  site  for  site, 
whereas  Event  5  has  a  delayed  phase,  possibly  a  pP  phase,  which  persists  across  all 
traces  of  the  array.  This  delayed  phase  is  not  obvious  at  LASA  nor  at  the  other  arrays. 

In  order  to  test  the  self-consistency  of  the  data,  transfer  functions  were 
first  computed  to  shape  Events  1,  2  and  4  to  Event  5  at  each  sensor  of  the  arrays. 

This  is  equivalent  to  computing  the  spectral  ratio  of  two  events  in  the  frequency 
domain,  including  the  phase  information.  Due  to  the  tight  cluster  of  epicenters,  the 
transmission  path  effects  from  each  source  to  a  given  receiver  are  assumed  to  be 

3 

equal,  thus  dividing  out  in  the  frequency  domain.  This  is  a  well  known  technique  for 
eliminating  unknown  but  common  transmission  effects  from  pairs  of  signals.  Such  a 
spectral  ratio  should  therefore  equal  the  spectral  ratio  of  the  source  radiations  for  the 
two  events  along  the  common  take  off  direction  from  source  to  receiver.  Unfortunately, 
each  source  radiation  includes  the  depth  of  burial  effect,  which  is  an  unknown  factor 
mixed  into  the  transfer  functions  of  the  data. 
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Fig.  6.  LASA  seismograms  of  the  four  events.  They  are  arranged  in  order 
of  increasing  LASA  magnitude  from  left  to  right.  The  relative  amplitudes 
from  subarray  to  subarray  of  each  event  are  preserved. 


SITE  NUMBER 


-2- 9772 


Fig.  7.  OONW  seismograms  of  three  events.  They  are  arranged 
in  order  of  increasing  amplitude  from  left  to  right.  The  relative 
amplitudes  from  sensor  to  sensor  of  each  event  are  preserved. 
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If  the  source  to  receiver  transmission  effects  are  indeed  common  for  events 
recorded  at  a  given  sensor,  then  sets  of  transfer  functions,  or  spectral  ratios,  should 
be  obtained  which  are  locally  consistent  at  LASA  and  at  OONW  in  spite  of  the  strong 
signal  variations  observed  from  sensor  to  sensor. 

Let  Elk(t)  be  the  seismogram  of  Event  i  recorded  at  sensor  k  of  either  array. 
In  the  frequency  domain  we  define  the  transfer  function  R1Jk(co)  as  the  ratio 


RijvcM 


E.ikM  ^  AM. 

Elk(o>)  S^co) 


(1) 


where  St  (to)  is  the  source  spectrum  of  Event  i  along  the  take  off  angle  to  the  array 
considered.  Transfer  functions  were  computed  in  time  by  a  least  squares  method  to 
shape  each  of  the  smaller  events  to  Event  5,  at  the  sites  of  each  array.  The  sampling 
increments  of  the  data  is  .05  sec.  In  each  case  a  50  point  transfer  function  R1Jk(t)  was 
computed  such  that  the  convolution  of  100  points  of  Ellc(t)  with  R1Jk(t)  is  the  best  least 
squares  approximation  to  149  points  of  EJk(t).  This  method  is  described  in  the  Appen¬ 
dix. 

Figures  8,  9  and  10  show  the  filters  Ri6ic(t),  R25k(t)  and  R46k(t)  respectively 
at  LASA.  Index  k  runs  from  1  to  21  corresponding  to  the  subarray  sum  traces  at  FI, 
F2,  .  .  .  El,  E2,  .  .  .  A0  respectively.  In  each  figure  the  average  filter  Rt }  (t)  is 
also  shown.  The  spike  at  the  end  of  some  filters  is  a  spurious  effect  which  should  be 
ignored  (see  Appendix). 

There  is  fairly  good  coherence  of  the  transfer  functions  R1Bk(t)  and  Es6k(t) 
from  sensor  to  sensor  despite  the  strong  signal  variations  across  LASA  seen  in  Figure 
8,  and  the  transfer  functions  R46k(t)  are  extremely  coherent. 
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Fig.  8.  Transfer  functions  Risk 00  at  each 
available  LAS  A  subarray.  The  average  filter 
Rls(t)  is  also  displayed. 


Fig.  9.  Transfer  functions  R-25K^  at  63 
available  LASA  subarray.  The  average  filter 
R25OO  is  also  shown. 
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Fig.  10.  Transfer  ^functions  at  ea°h  available  LASA  subarray. 

The  average  filter  £45(1)  is  also  shown. 
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In  Figures  11  and  12  the  transfer  functions  Ris^t)  and  R^ft)  are  displayed 
for  the  available  data  at  OONW.  These  transfer  functions  are  also  very  coherent  but 
are  more  oscillatory  than  the  LASA  transfer  functions. 

Similar  calculations  using  the  individual  sensors  of  the  United  Kingdom 
arrays  were  done.  Locally  consistent  transfer  functions  were  also  obtained  at  these 
arrays,  although  the  signal  to  noise  ratio  of  the  data  was  not  as  good  as  the  LASA  data. 
From  this  one  can  conclude  that  variations  in  the  source  radiation  of  these  events  are 
not  detected  across  the  aperture  of  each  array  and  that  transmission  path  effects  do 
cancel  out.  This  implies  that  steered  beams  of  each  event  can  be  used  to  calculate 
transfer  functions  without  degrading  source  information. 
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Fig.  11.  Transfer  functions  R,c^(t) 
at  each  available  sensor  at  OOnW. 
R^5(t)  is  the  average  filter. 
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Fig.  12.  Transfer  functions  Rosk^) 
at  each  available  sensor  at  OONW. 
R25(t)  is  the  average  filter. 


IV.  TRANSFER  FUNCTIONS  USING  BEAM  DATA 


Figure  13  shows  the  steered  beam  traces  at  each  array  except  LASA.  At 

LASA  the  straight  sum  trace  of  the  F3  subarray  is  shown  in  order  to  compare  arrays 

of  roughly  the  same  size.  Although  there  is  some  degradation  of  high  frequencies  in 

4 

the  straight  sum  trace  ,  which  could  be  removed  by  beam  steering  the  individual 
traces,  this  degradation  is  common  to  all  four  events  and  will  not  affect  the  data 
analysis.  Event  5  is  normalized  in  Figure  6  to  have  the  same  trace  amplitude  at  each 
array.  The  other  three  events  are  plotted  with  the  correct  amplitudes  relative  to 
Event  5.  The  local  array  magnitude  for  each  event  is  indicated,  those  magnitudes 
with  asterisks  (*)  being  estimated  by  comparing  amplitudes  with  other  events  at  the 
same  array.  The  instrument  calibrations  at  the  OONW  array  are  unfortunately  not 
known,  but  the  relative  amplitudes  shown  for  the  three  events  are  considered  reliable. 

A  glance  at  these  beams  shows  that  the  seismic  radiation  patterns  of  the  four 
events  are  not  equal.  At  LASA,  for  example,  there  is  a  magnitude  increase  of  .  7 
units  from  Event  1  to  Event  5,  whereas  at  GBA,  Events  1  and  5  have  the  same  magnitude. 
Another  strange  feature  is  the  apparent  discrepancy  in  radiation  patterns  recorded  at 
LASA  and  YKA,  although  the  take  off  angles  and  azimuths  from  the  source  region  are 
nearly  identical  to  both  arrays.  The  precursors  seen  at  YKA  on  Events  1  and  2,  which 
are  clearly  seen  on  the  individual  sensors,  are  not  present  on  Events  4  and  5.  At 
LASA,  no  such  precursors  are  seen  on  Events  1  and  2.  At  OONW  an  apparent  depth 
phase  occurs  on  Event  5  which  is  not  seen  on  Events  1  and  2.  Yet  at  GBA,  which  is  the 
same  distance  from  the  Soviet  Test  Site  but  along  a  different  azimuth.  Events  5  and  1 
are  quite  similar  with  no  apparent  depth  phase  seen. 
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Fig.  13.  Array  beams  of  the  four  events.  At  LASA  the  straight  sum  for 
subarray  F3  is  shown.  All  other  beams  are  steered. 
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These  discrepancies  in  amplitude  and  wave  shape  are  probably  caused  by 
site  conditions  which  are  difficult  to  separate.  Small  changes  in  epicenter  and  burial 
depth,  the  orientation  of  the  initial  stress  field  in  the  source  region,  and  non -sphericity 
of  the  explosive  pressure  pulse  are  likely  factors  for  producing  anisotropic  radiation 
to  teleseismic  distances. 

Transfer  functions  Ris(t),  R25  (t)  and  R46(t)  were  computed  at  each  array 
using  the  beam  data  of  Figure  13.  These  functions  are  displayed  in  Figure  14.  If  the 
seismic  source  radiations  for  each  event  were  indeed  isotropic,  then  the  functions  down 
each  column  would  have  the  same  shape  and  amplitude  since  transmission  path  effects 
have  been  eliminated.  R45(t)  seems  to  be  consistent  from  array  to  array,  whereas 
R2  5  (t)  and  R15(t)  increase  in  amplitude  with  epicentral  distance  from  GBA  to  LASA, 
and  then  an  amplitude  drop  to  WRA.  Thus  Events  4  and  5  seem  to  have  similar  radia¬ 
tion  patterns,  however  complex,  but  Events  1  and  5  clearly  do  not.  The  spikes  at  the 
end  of  some  transfer  functions  should  be  ignored  (see  Appendix). 

Of  the  five  arrays  YKA,  LASA  and  WRA  each  have  a  complete  set  of  three 
transfer  functions.  The  numerical  quality  of  the  transfer  functions  is  given  by  the 
normalized  error  e^  described  in  the  Appendix,  ejjj  =  0  indicates  a  perfect  filter  with 
no  error,  and  e^  =  1  is  the  worst  case  in  which  the  least  squares  filter  is  zero.  To 
illustrate  this  quality.  Event  1  was  convolved  with  Ris(t)  and  compared  with  Event  5  at 
arrays  YKA,  LASA  and  WRA.  These  results  are  shown  in  Figure  15.  Clearly,  the 
convolution  output  for  the  YKA  data  is  a  poor  approximation  to  Event  5,  whereas  at 
LASA  the  filtered  data  is  in  excellent  agreement  with  Event  5.  The  normalized  error 
ejjj  is  .  38  at  YKA,  .  06  at  LASA  and  .  17  at  WRA.  In  Figure  14  all  transfer  functions, 
except  those  at  YKA,  perform  as  well  or  better  than  the  WRA  example  shown  in 
Figure  15. 
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Fig.  14.  Transfer  functions  Rj^(t)  at  each  available  array  for  i  =  1,  2,  4. 
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ll»-2-9ng| 


Y  K  A  LASA  WRA 


Fig.  15.  Examples  of  numerical  accuracy  of  computed  transfer  functions  Rj^t). 
In  each  column  the  middle  trace  is  the  convolution  of  Event  1  with  Rj^t),  which 
can  be  compared  with  Event  5,  the  bottom  trace.  Rj^t)  at  LASA  is  quite  good, 
but  the  filter  at  YKA  is  poor  quality. 
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V.  THEORETICAL  TRANSFER  FUNCTIONS 


In  spite  of  the  apparent  anisotropy  of  the  source  radiations  for  the  events 


studied,  it  is  of  interest  to  calculate  theoretical  transfer  functions  for  assumed 

explosion  models  and  see  if  such  calculations  match  any  trends  in  the  observed  LASA 

transfer  functions,  which  were  the  most  reliable,  numerically. 

Two  isotropic  mathematical  models  for  the  shape  of  compressional  waves 

radiated  from  explosive  sources  are  those  of  Blake"’  and  Haskell^.  Blake  generalized 

7 

a  solution  obtained  by  Sharpe  for  P  waves  radiating  from  a  cavity  in  an  infinite  homo¬ 
geneous  elastic  medium,  where  the  cavity  is  excited  by  a  step  function  of  pressure. 
Haskell  took  the  near  in  displacement  records  of  nuclear  explosions  from  different 

g 

media  obtained  by  Werth  and  Herbst  and  fit  the  record  shape  by  a  family  of  analytic 
functions  which  could  be  adjusted  for  different  media  and  explosive  yield. 

Blake's  expression  for  the  far  field  radial  particle  velocity,  S^t),  produced 
by  applying  a  step  function  of  pressure  pt  to  an  elastic  cavity  wall  is 


(2) 


where 


r 


radial  distance  from  source  to  receiver 


c 


compressional  velocity  of  medium 


CT 


Poisson's  ratio 


p  =  density  of  medium 
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a,  =  (c/a,)  a  -  2a)2/ a  -  CT) 

w,  =  (c/a,)  (1  -  2ct)/(1  -  a)  =  (1  -  2a)2a, 


and 


-i  i 

0  =  tan  (1  -  2g)2 


For  explosions  the  radius  a  is  taken  to  be  the  equivalent  elastic  radius, 
outside  of  which  the  medium  behaves  elastically,  whereas  pressure  p  is  determined  by 
the  lithostatic  pressure  at  the  test  site.  Both  of  these  variables  are  functions  of  the 
burial  depth  which  is  not  known  for  these  events.  As  noted  by  Filson\  cepstral 
analysis  of  these  four  events  does  not  reveal  consistent  delay  times  at  the  five  arrays 
for  possible  pP  phases.  We  shall  assume  that  all  four  events  have  equal  depths  in 
which  case  the  pressures  p,  assumed  for  each  elastic  cavity  are  equal.  For  the  com¬ 
putation  of  theoretical  transfer  functions  for  events  recorded  at  the  same  site,  it  is 
sufficient  to  write  the  source  function  as 


S,(t)  ~  a,e  ait  cos  (oo  ,t  +  0) 


(3) 


where  the  tilda  (~)  denotes  proportionality. 

Another  source  model  we  shall  consider  is  obtained  from  Haskell's  reduced 
displacement  potential  i|r .  In  this  case  the  far  field  radial  particle  velocity  is  given  by 


■  TF  ^ 


which  for  events  recorded  at  the  same  site  can  be  written 
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/ 


s^t)  ~  ti  Hk^e"klt  (1  +  24B)  - 


(^3 


0  +  48B)  +  B(ktt)4 


(4) 


The  parameters  f  (®),  k  and  B  were  used  by  Haskell  to  generate  a  family  of  analytic 
functions  t|r (t)  which  would  scale  properly  with  yield  Y  (kilotons).  For  scaling  purposes 
Haskell  assumed  that 


♦  («)  ~  Y  (5) 

k  ~  Y~3  (6) 

and 

a  =  100  Y^  (7) 

where  a  is  the  equivalent  elastic  radius  in  meters,  from  the  source,  outside  of  which 
the  medium  behaves  elastically.  This  radius  will  be  used  in  equation  (3)  for  Blake's 
far  field  solution. 

Under  these  gross  assumptions  the  source  function  for  either  model,  equation 

1 

(3)  or  (4),  has  amplitude  and  time  scales  which  are  proportional  to  Y3  at  a  given  array. 

Yields  are  not  available  for  these  events;  however,  an  approximate  magnitude  yield 

relationship  can  be  assumed.  The  Soviet  Test  Site  is  within  the  northeastern  edge  of 

the  Balkhash  Chingiz  Foldbelt,  a  region  of  intrusive  igneous  rocks  overlain  by  folded 

sediments.  We  shall  assume  that  the  source  medium  for  each  shot  is  granite.  For 

1  9 

explosions  in  hard  rock  there  is  evidence  ’  that  body  wave  magnitude  up  to  about 
6.  can  be  approximated  by  the  formula 
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(8) 


mb  =  3.  8  +  log1Q  Y 

This  differs  slightly  from  other  magnitude  yield  relationships,  e.g.  those 
summarized  by  Marshall^,  but  these  differences  are  negligible  for  this  analysis, 
considering  the  lack  of  consistency  from  array  to  array. 

From  the  LASA  body  wave  magnitudes  a  set  of  yields  were  calculated  for 
the  four  events  using  equation  (8).  These  are  listed  in  Table  2.  From  these  yields, 
values  of  i|f  (=°)  and  k  were  computed  from  (5)  and  (6)  by  scaling  up  Haskell’s  values  of 
k  =  31.  6  sec  1  and  i|t  («>)  =  2.  5  x  103m3  for  a  5  kt  explosion  in  granite.  These  parameters 
and  values  of  at  calculated  from  (7)  are  also  tabulated. 

Table  2 


Event 

LASA  mb 

Y, 

k, 

Ji-H 

at 

1 

5.4 

39.8  kt 

15.  8  sec  1 

19.  4  x  lO^3 

341  m 

2 

5.6 

63.0 

13.6 

31.  5  x  103 

398 

4 

5.8 

100. 

11.  7 

50.  0  x  103 

464 

5 

6.1 

199. 

9.16 

100.  x  103 

584 

Yield  dependent  parameters  used  in  source  functions  for  the  events  recorded  at  LASA. 

The  remaining  parameters  required  for  calculating  the  source  functions  are 
independent  of  yield,  and  depend  on  the  source  medium  only.  For  Blake's  far  field 
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solution  values  of  c  =  5  km/sec,  and  t=  .  3  were  assumed,  and  in  Haskell's  formula  a 
value  of  the  parameter  B  =  .  24  for  granite  was  used. 

From  the  above  numbers  source  functions  were  calculated  for  each  event 
using  equations  (3)  and  (4).  These  functions  are  shown  in  Figure  16.  The  main  differ¬ 
ence  between  the  two  source  functions  is  that  a  step  function  of  pressure  in  Blake's 
solution  causes  a  discontinuity  in  the  velocity  response  at  zero  time,  whereas  Haskell's 
analytical  function  was  constrained  to  have  continuous  displacement,  velocity,  and 
acceleration  at  zero  time. 

Theoretical  transfer  functions  Ris(t),  R26(t)  and  R4B(t)  were  calculated  for 
both  sets  of  assumed  source  functions.  These  filters  are  superposed  on  the  observed 
LASA  transfer  functions  in  Figure  17. 


|  H-2-95SH-I  | 


Fig.  16.  Particle  velocity  source  functions  for  the  four  events 
using  Blake’s  and  Haskell's  model. 
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LASA  ( F 3)  -  OBSERVED 

BLAKE  - 

HASKELL 


} 


THEORETICAL 


Fig.  17.  Comparison  of  theoretical  and  observed  transfer  functions  at  LASA. 


28 


VI.  DISCUSSION  OF  RESULTS  AND  CONCLUSIONS 


From  Figure  17  it  appears  that  the  transfer  functions  obtained  from  Blake’s 
solution  fit  the  observed  LASA  data  better  than  those  obtained  from  Haskell’s  source 
functions.  Disregarding  the  2-3  Hz  oscillations,  the  first  negative  swing  seen  in 
Ri  6  (t)  diminishes  in  R26(t)  and  R4s(t)  in  both  the  data  and  the  theoretical  filters  for 
Blake's  source  model. 

Both  source  models  predict  a  relative  degradation  of  high  frequencies  with 
increasing  magnitude.  This  implies  that  Rib  (to)  should  show  a  more  rapid  attenuation 
of  high  frequencies  than  either  R26(to)  or  R4B(to).  This  effect  is  demonstrated  in 
Figure  18,  which  shows  the  amplitude  spectra  of  the  three  LASA  transfer  functions. 

The  curves  have  been  normalized  to  unity  at  1  Hz  to  emphasize  the  different  attenuation 
rates  of  the  filters. 

A  major  drawback  in  this  analysis  is  that  the  depth  of  focus  of  each  event  has 
been  ignored.  Unfortunately,  the  depth  effect  might  contaminate  that  part  of  the  trans¬ 
fer  function  which  was  interpreted  in  terms  of  scaled  source  functions.  Taking  the 
simplest  example,  let  us  assume  that  Event  1  was  at  zero  depth,  and  Event  5  had  a 
depth  corresponding  to  a  pP-P  time  of  r  seconds.  Then  for  identical  point  sources,  the 
observed  transfer  function  at  teleseismic  distance  would  be  approximately 

Rie(t)  ~  6(t)  +  r6(t  -  t) 

where  r  is  the  plane  wave  reflection  coefficient  at  the  free  surface  and  6(t)  is  the  Dirac 
delta  function.  For  any  reasonable  crustal  model  r  varies  from  about  -.  7  to  - 1.  0 
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RELATIVE  AMPLITUDE 


Fig.  18.  Amplitude  response  of  the  LASA  transfer  functions. 
Curves  are  normalized  to  have  amplitude  1  at  1  Hz. 
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depending  on  the  take  off  angle  from  source  to  array.  If  Tis  taken  to  be  .  5  sec. ,  then 
RibW  will  consist  of  a  positive  spike  followed  by  a  smaller  negative  spike  .  5  sec. 
later.  This  operator,  if  band  pass  filtered,  would  have  a  shape  similar  to  the  first 
second  of  the  LASA  transfer  function  Rie(t)  shown  in  Figure  17. 

The  radiation  of  pP  to  teleseismic  distances  is  obviously  more  complicated, 
as  the  data  indicates.  Numerical  studies  by  Plamondon^  show  that  strong  multiple 
reflections  occur  between  the  free  surface  and  a  buried  spherical  cavity  which  has 
been  excited  by  a  pressure  pulse.  From  this  it  follows  that  the  observed  teleseismic 
radiation  from  shallow  cavities  cannot  be  simply  modelled  by  point  sources.  Clearly, 
numerical  studies  similar  to  those  of  Plamondon  are  needed  to  calculate  the  seismic 
radiation  along  different  take  off  angles  from  well  documented  explosions.  This  may 
explain  the  observed  anisotropy  in  source  radiation  and  enable  one  to  separate  the 
effects  of  burial  depth  and  seismic  scaling. 
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APPENDIX 


The  method  of  least  squares  filters  is  very  well  known.  One  of  the  earliest  and 

most  readable  references  for  practical  computations  using  discrete  time  data  is  by 
12 

Levinson  .  This  reference  also  appears  as  Appendix  B  of  N.  Wiener's  book  Time 
Series  ,  M.  I.  T.  Press,  1964. 

One  of  the  problems  in  the  calculations  of  this  paper  is  a  spike  which  sporadically 
occurs  at  the  ends  of  the  transfer  functions  (or  filters).  This  appears  to  be  a  function 
of  the  finite  length  data  windows  used  in  the  filter  calcuations.  In  order  to  discuss  this 
it  is  necessary  to  examine  the  matrix  equations  used  in  the  computations.  To  simplify 
the  discussion  let  us  consider  the  following  small  system  of  equations: 


yi 

0 

0 

y2 

x2 

xx 

0 

ys 

= 

x3 

x2 

Xi 

y4 

x4 

x3 

x2 

yS 

0 

x4 

X3 

Ye 

0 

m 

0 

X4 

This  equation  gives  the  output,  y  ,  obtained  by  convolving  the  input  data,  x^, 
with  a  filter  f^..  For  such  a  transient  convolution  the  number  of  samples  in  y  is  one 
less  than  the  sum  of  the  filter  and  data  lengths,  i. e.  3  +  4  -  1  =  6.  Thus,  if  y  and  x^. 
are  given,  the  system  of  equations  is  overdetermined  for  calculating  f  .  However  a 
solution  for  f  can  be  obtained  which  minimizes  the  sum  of  squared  errors  between  the 
desired  output,  y  ,  and  the  actual  convolution. 

If  we  write  (Al)  in  the  following  matrix  notation 
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y 


=  Xf  ,  (A2) 

then  premultiplying  (Al)  by  XT  gives 


XTXf  = 

XTy 

(A3) 

ro 

ri 

r3 

where 

XTX  = 

rl 

To 

ri 

=  R, 

r3 

rl 

r0 

and 

4-i+l 

rT  - 
T  i=l 

»  T 

=  0, 

1,  2. 

(A4) 

The  solution  for  f  in  (A3)  is  the  least  squares  solution  of  (Al).  R  is  a  positive  definite 

matrix  and  each  element  r^  is  the  transient  autocorrelation  of  the  input  data  at  lag 

t.  The  elements  along  any  diagonal  of  Rare  equal,  hence  R  is  called  a  Toeplitz  matrix. 

13 

Large  matrices  of  this  form  can  be  inverted  very  rapidly  using  Levinson’s  algorithm  . 

Referring  to  Section  IV  of  this  paper  we  let  xt  be  Ex(t),  and  y  be  Es(t)  so  that 
ft  corresponds  to  the  transfer  function  Ris(t)  computed  at  one  of  the  arrays.  In  order 
to  apply  an  equation  like  (Al)  to  real  data,  E^t)  and  E3(t)  must  be  truncated  even  though 
it  is  often  not  clear  when  a  P  wave  terminates.  From  equation  (Al)  one  sees  that  ys 
and  y6  are  not  functions  of  X5  and  x6  because  only  four  points  of  x  were  used.  This 
may  introduce  numerical  errors  into  f3  and  f3  which  are  physically  not  meaningful. 
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One  possible  test  of  the  sensitivity  of  f2  and  f3  to  this  data  is  to  set  y6  =  Ye  =  0  in  (Al). 
This  will  have  no  effect  on  fx ,  but  may  change  f3  and  f3  significantly. 

This  idea  was  checked  for  a  large  system  of  equations  similar  to  (Al)  using  Ex  (t) 
and  E5(t)  recorded  at  LASA.  In  Figure  19  results  of  a  filter  calculation  are  shown 
using  60  samples  of  E^t)  as  input,  and  109  points  of  Es(t)  as  desired  output.  A  least 
squares  transfer  function  50  samples  long  was  computed.  This  filter  shows  a  large 
positive  spike  at  the  end. 

In  Figure  20  the  same  length  filter  was  computed,  but  in  this  case  all  samples  of 
Es(t)  past  number  60  were  arbitrarily  set  equal  to  zero.  In  this  case  the  spike  at  the 
end  of  the  filter  has  been  removed,  but  the  early  part  of  the  filter  remains  unchanged. 
The  convolution  of  the  filter  with  the  input  data  yields  a  good  approximation  to  the 
desired  output,  even  for  the  zeroed  portion  of  the  data. 

The  error  minimized  in  designing  the  filter  is  the  sum  of  squares 


e' 


.2 


(y  -  Xf  )T  (y  -  Xf  ) 


(A5) 


Using  (A3),  this  reduces  to 


yTTy  -  f  TXX  f 


(A6) 


Dividing  e2  by  the  sum  of  squares  yTy  we  obtain  the  normalized  error 


e2  =  1  -  f T  XXf 


(A  7) 
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Fig.  19.  Transfer  function  Rjc(t)  at  LASA 
computed  using  60  points  of  Ei(t)  and  109 
points  of  Eij(t).  The  convolution  Ej^t^Rj^t) 
is  a  good  approximation  to  Eg(t). 


Fig.  20.  Transfer  function  Ris(t)  computed 
as  in  Fig.  19  but  with  only  60  points  of  E5(t) 
being  used.  R^t)  shows  no  spike  at  the 
end  in  this  case. 


which  falls  in  the  range  0  to  1.  If  ejj  =  0,  the  filter  exactly  converts  xt  to  y  ,  and  if 
e^  =  1,  the  filter  does  the  poorest  job  of  shaping  x^  to  y  .  This  occurs  when  x^  and  y{ 
are  uncorrelated  and  all  f  are  zero. 

The  filters  computed  at  LASA,  for  example,  all  have  normalized  errors  of 
<0. 1.  At  the  other  arrays  the  filters  have  normalized  errors  as  large  as  0.  3. 

The  analysis  in  this  paper  of  the  LASA  filters  in  terms  of  source  functions  was 
pertinent  only  to  the  first  20  points  (1  second)  of  the  filters.  These  filter  points  are 
numerically  stable  and  show  very  little  dependence  on  the  length  of  x  ,  y  or  f^ 
specified  in  the  calculations,  as  long  as  f^  has  more  than  about  30  points. 

The  effect  of  the  convolution  "tail"  in  (Al)  can  be  eliminated  by  using  the  same 
number  of  samples  in  x^  and  y  .  This  can  be  done  by  taking  only  the  first  four  equations 
of  (Al),  for  example,  and  repeating  the  analysis  through  (A7).  The  only  difficulty  in 
this  case  is  that  R  is  no  longer  a  Toeplitz  matrix  and  is  considerably  more  time- 
consuming  to  invert. 
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13.  ABSTRACT 


Short  period  P  waves  of  four  presumed  Soviet  explosions  from  Eastern  Kazakh  are  examined 
at  five  teleseismic  arrays:  LASA,  YKA,  OONW,  WRA  and  GBA.  Transfer  functions  to  shape  the 
lower  magnitude  to  the  highest  magnitude  event  were  computed  at  each  array  to  eliminate  trans¬ 
mission  path  effects  from  source  to  receiver.  The  transfer  functions  are  locally  consistent  from 
sensor  to  sensor  at  each  array,  but  show  considerable  global  variations  from  array  to  array. 

This  suggests  that  there  are  strong  azimuthal  variations  in  the  source  radiation  of  the  events,  due 
to  complex  scattering  of  the  signals  at  the  test  site.  At  LASA,  the  observed  transfer  functions 
can  be  explained  by  using  explosion  source  functions  of  Blake  and  Haskell.  The  assumed  source 
parameters  are  scaled  functions  of  the  yield  of  each  event,  which  is  estimated  from  an  empirical 
amplitude  yield  relationship. 
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